home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / matrix.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  22.2 KB  |  688 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module matrix)
  13.  
  14. (DECLARE-TOP(SPECIAL ERRRJFFLAG ONEOFF* EI* EJ* *ECH* *TRI* *INV*
  15.           MDL DOSIMP $DETOUT VLIST MUL* TOP* *DET* GENVAR $RATFAC
  16.           *MOSESFLAG VARLIST HEADER LININD* $SCALARMATRIXP $SPARSE
  17.           $ALGEBRAIC *RANK*) 
  18.      (*LEXPR FMAPL1) (FIXNUM NN LEN)
  19.      (GENPREFIX X))
  20.  
  21. (DEFMVAR $DETOUT NIL)
  22. (DEFMVAR TOP* NIL)
  23. (DEFMVAR $RATMX NIL)
  24. (DEFMVAR $MATRIX_ELEMENT_MULT '|&*|)    ;;; Else, most useful when '|&.|
  25. (DEFMVAR $MATRIX_ELEMENT_ADD '|&+|)
  26. (DEFMVAR $MATRIX_ELEMENT_TRANSPOSE NIL)
  27.  
  28. ;;provides some of the same spirit of *array but
  29. ;;in the value cell. see get-array-pointer below.
  30. (defun cl-*array (nam maclisp-type &rest dimlist)
  31.   (proclaim (list 'special nam))
  32.   (set nam (apply '*array nil maclisp-type  dimlist)))
  33.  
  34. #+MacLisp
  35. (DEFUN GET-ARRAY-POINTER (X)
  36.   (COND ((EQ (ml-typep X) 'array) X)
  37.     ((GET X 'array))
  38.     (T (MERROR "~S is not an array." X))))
  39.  
  40. #+Franz
  41. (defun get-array-pointer (x)
  42.    (cond ((arrayp x) x)
  43.      ((and (symbolp x) (arrayp (getd x))) x)
  44.      (t (merror "~s is not an array." x))))
  45.  
  46. #+oldlispm
  47. (DEFUN GET-ARRAY-POINTER (X)
  48.   (COND ((ARRAYP X) X)
  49.     ((FBOUNDP X) (SYMBOL-FUNCTION X))
  50.     (T (ERROR  "~S is not an array." X))))
  51. #+NIL    ;Defined by maclisp-compatibility array stuff, for just this purpose.
  52. (deff get-array-pointer
  53.   #'si:get-array-pointer)
  54.  
  55. ;;I believe that all the code now stores arrays in the value cell 
  56. (DEFUN GET-ARRAY-POINTER (SYMBOL)
  57.   "There may be nesting of functions and we may well need to apply
  58.    this twice in a row"
  59.   (cond ((arrayp symbol) symbol) (t (symbol-value symbol))))
  60.  
  61. (DEFUN MXC (X) (MAPCAR #'(LAMBDA (Y) (CONS '(MLIST) Y)) X))
  62.     ; Matrix to MACSYMA conversion
  63.  
  64. (DEFUN MCX (X) (MAPCAR #'CDR X))  ; MACSYMA to Matrix conversion
  65.  
  66. (DEFUN TRANSPOSE (M)
  67.        (PROG (B NN LEN)
  68.          (SETQ LEN (LENGTH (CAR M)) NN 1)
  69.     LOOP (COND ((> NN LEN) (RETURN B))) 
  70.          (SETQ B (NCONC B (NCONS (NTHCOL M NN))) NN (f1+ NN))
  71.          (GO LOOP)))
  72.  
  73. (DEFUN NTHCOL (X NN)
  74.        (COND ((OR (NULL X) (> NN (LENGTH (CAR X)))) NIL) (T (NTHCOL1 X NN))))
  75.  
  76. (DEFUN NTHCOL1 (X NN)
  77.        (COND ((OR (NULL X) (= NN 0)) NIL)
  78.          (T (CONS (ITH (CAR X) NN) (NTHCOL1 (CDR X) NN)))))
  79.  
  80. (DEFUN CHECK (X) (COND ((ATOM X) (MERROR "Not matrix:~%~M" X))
  81.                ((EQ (CAAR X) '$MATRIX) X)
  82.                ((EQ (CAAR X) 'MLIST) (LIST '($MATRIX) X))
  83.                (T (MERROR "Not matrix:~%~M" X)))) 
  84.  
  85. (DEFUN CHECK1 (X) (COND ((ATOM X) NIL)
  86.             ((EQ (CAAR X) '$MATRIX) X)
  87.             ((EQ (CAAR X) 'MLIST) (LIST '($MATRIX) X)))) 
  88.  
  89. (DEFMFUN $MATRIXP (X) (AND (NOT (ATOM X)) (EQ (CAAR X) '$MATRIX)))
  90.  
  91. (DEFMFUN $CHARPOLY (MAT VAR) 
  92.        (SETQ MAT (CHECK MAT))
  93.        (IF (NOT (= (LENGTH MAT) (LENGTH (CADR MAT))))
  94.        (MERROR "Matrix must be square - CHARPOLY")) 
  95.        (COND ((NOT $RATMX) (DET1 (ADDMATRIX1
  96.                 (SETQ MAT (MCX (CDR MAT))) 
  97.                 (DIAGMATRIX (LENGTH MAT) (LIST '(MTIMES) -1 VAR) '$CHARPOLY))))
  98.          (T (NEWVAR VAR) (NEWVARMAT1 MAT)
  99.         (SETQ MAT (MCX (CDR MAT)))
  100.         (DETERMINANT1 (ADDMATRIX MAT (DIAGMATRIX (LENGTH MAT) 
  101.                              (LIST '(MTIMES) -1 VAR)
  102.                              '$CHARPOLY))))))
  103.  
  104. (DEFUN DISREPLIST1 (A) (SETQ HEADER (LIST 'MRAT 'SIMP VARLIST GENVAR))
  105.                (MAPCAR #'DISREPLIST A))
  106.  
  107. (DEFUN DISREPLIST (A) (MAPCAR #'(LAMBDA (E) (CONS HEADER E)) A))
  108.  
  109. (DEFUN REPLIST1 (A) (MAPCAR #'REPLIST A)) 
  110.  
  111. (DEFUN REPLIST (A) (MAPCAR #'(LAMBDA (E) (CDR (RATREP* E))) A))
  112.  
  113.  
  114. (DEFUN TIMEX (MAT1 MAT2)
  115.  (COND ((EQUAL MAT1 1) MAT2)
  116.        ((AND ($MATRIXP MAT1) ($MATRIXP MAT2) (NULL (CDR MAT1)))
  117.     (NCONS '($MATRIX SIMP)))
  118.        (T (NEWVARMAT MAT1 MAT2)
  119.       (LET (($SCALARMATRIXP
  120.          (IF (AND ($LISTP MAT1) ($LISTP MAT2)) T $SCALARMATRIXP)))
  121.            (SIMPLIFYA (TIMEX0 MAT1 MAT2) NIL)))))
  122.  
  123. (DEFUN LNEWVAR (A)
  124.        ((LAMBDA (VLIST)
  125.         (LNEWVAR1 A)
  126.         (SETQ VARLIST (NCONC (SORTGREAT VLIST) VARLIST)))
  127.     NIL))
  128.  
  129. (DEFUN LNEWVAR1 (A)
  130.        (COND ((ATOM A) (NEWVAR1 A))
  131.          ((MEMQ (CAAR A) '(MLIST MEQUAL $MATRIX)) (MAPC #'LNEWVAR1 (CDR A)))
  132.          (T (NEWVAR1 A))))
  133.  
  134. (DEFUN NEWVARMAT (MAT1 MAT2)
  135.        (COND ($RATMX
  136.           ((LAMBDA (VLIST)
  137.         (LNEWVAR1 MAT1) (LNEWVAR1 MAT2)
  138.         (SETQ VARLIST (NCONC (SORTGREAT VLIST) VARLIST))) NIL))))
  139.  
  140. (DEFUN NEWVARMAT1 (A)
  141.        (COND ($RATMX (LNEWVAR A))))
  142.  
  143. (DEFUN ADDMATRIX (X Y) (SETQ X (REPLIST1 X) Y (REPLIST1 Y))
  144.                (DISREPLIST1 (ADDMATRIX1 X Y)))
  145.  
  146. (DEFUN ADDMATRIX1 (B C)
  147.        (COND ((NOT (AND (= (LENGTH B) (LENGTH C))
  148.             (= (LENGTH (CAR B)) (LENGTH (CAR C)))))
  149.           (MERROR "Attempt to add stuff of unequal length")))
  150.        (MAPCAR #'ADDROWS B C))
  151.  
  152. (DEFUN ADDROWS (A B)
  153.        (COND ((NOT $RATMX) (MAPCAR #'(LAMBDA (I J)
  154.                          (SIMPLUS (LIST '(MPLUS) I J) 1 NIL)) A B))
  155.          (T (MAPCAR #'RATPLUS A B)))) 
  156.  
  157. (DEFMFUN $DETERMINANT (MAT)
  158.   (COND ((ATOM MAT) (LIST '(%DETERMINANT) MAT))
  159.     (T (SETQ MAT (CHECK MAT))
  160.        (IF (NOT (= (LENGTH MAT) (LENGTH (CADR MAT))))
  161.            (MERROR "DETERMINANT called on a non-square matrix."))
  162.            (COND ((NOT $RATMX) (DET1 (MCX (CDR MAT))))
  163.              (T (NEWVARMAT1 MAT) (DETERMINANT1 (MCX (CDR MAT))))))))
  164.  
  165. (DEFUN DET (M)
  166.   (IF (= (LENGTH M) 1)
  167.       (CAAR M)
  168.       (LET (*DET* MUL*)
  169.     (MTOA '*MAT* (SETQ *DET* (LENGTH M)) *DET* M)
  170.     (SETQ *DET* (TFGELI0 '*MAT* *DET* *DET*))
  171.     (RATREDUCE *DET* MUL*)))) 
  172.  
  173. (DEFUN DETERMINANT1 (X) (CATCH 'DZ (RDIS (DET (REPLIST1 X))))) 
  174.  
  175. (DEFUN TREEDET (MAT)
  176.   (PROG (ROW MDL LINDEX TUPLEL N ID MD LT)
  177.     (SETQ MAT (REVERSE MAT))
  178.     (SETQ N (LENGTH MAT) MD (CAR MAT))
  179.     (SETQ MAT (CDR MAT))(SETQ LINDEX (NREVERSE (INDEX* N)) TUPLEL (MAPCAR #'LIST LINDEX))
  180.      LOOP1(COND ((NULL MAT) (RETURN (CAR MD))))
  181.     (SETQ MDL NIL)
  182.     (MAPCAR #'(LAMBDA(A B)
  183.              (SETQ MDL(NCONC MDL (LIST A B))))
  184.         TUPLEL MD)
  185.     (SETQ MD NIL)
  186.     (SETQ ROW (CAR MAT)MAT (CDR MAT))
  187.     (SETQ LT (SETQ TUPLEL (NEXTLEVEL TUPLEL LINDEX)))
  188.      LOOP2(COND ((NULL LT) (SETQ MD (NREVERSE MD)) (GO LOOP1)))
  189.     (SETQ ID (CAR LT) LT (CDR LT)) (SETQ MD (CONS (COMPUMD ID ROW) MD)) (GO LOOP2) ))
  190.  
  191. (DEFUN ASSOO (E L) (PROG()
  192.              LOOP(COND ((NULL L) (RETURN NIL))
  193.                    ((EQUAL E (CAR L)) (RETURN (CADR L))))
  194.             (SETQ L (CDDR L))(GO LOOP)))
  195.  
  196. (DEFUN COMPUMD (ID ROW)
  197.   (PROG(E MINOR I D SIGN ANS)
  198.        (SETQ ANS 0 SIGN -1 I ID)
  199.     LOOP(COND ((NULL I)(RETURN ANS))) 
  200.        (SETQ D (CAR I) I (CDR I) SIGN (TIMES -1 SIGN))
  201.        (COND ((EQUAL (SETQ E(ITH ROW D)) 0)(GO LOOP))
  202.          ((EQUAL (SETQ MINOR(ASSOO (zl-DELETE D(COPY ID)) MDL)) 0)(GO LOOP)))
  203.        (SETQ ANS (SIMPLUS (LIST '(MPLUS) ANS (SIMPTIMES (LIST '(MTIMES) SIGN E MINOR) 1 NIL)) 1 NIL)) (GO LOOP)))
  204.  
  205. ;Gag me with a vax!  --gsb
  206. ;(DECLARE(SPECIAL LTP*))
  207. ;
  208. ;(DEFUN APDL (L1 L2)
  209. ;  ((LAMBDA (LTP*)
  210. ;     (MAPCAR #'(LAMBDA (J) (SETQ LTP* (CONS (APPEND L1 (LIST J)) LTP*))) L2)
  211. ;     (NREVERSE LTP*))
  212. ;   NIL))
  213. ;
  214. ;(DECLARE(UNSPECIAL LTP*))
  215.  
  216. (defun apdl (l1 l2)
  217.   (mapcar #'(lambda (j) (append l1 (list j))) l2))
  218.  
  219. (DEFUN NEXTLEVEL (TUPLEL LINDEX)
  220.   (PROG(ANS L LI)
  221.     LOOP (COND ((NULL TUPLEL )(RETURN ANS)))
  222.        (SETQ L (CAR TUPLEL) TUPLEL (CDR TUPLEL) LI (CDR (NCDR LINDEX (CAR (LAST L)))))
  223.        (COND ((NULL LI) (GO LOOP)))
  224.        (SETQ ANS(NCONC ANS (APDL L LI))) (GO LOOP)))
  225.  
  226. (DEFUN DET1 (X)
  227.   (COND ($SPARSE (MTOA '*MAT* (LENGTH X) (LENGTH X) 
  228.                (MAPCAR #'(LAMBDA (X) (MAPCAR #'(LAMBDA (Y) (NCONS Y)) X))X))
  229.          (SPRDET '*MAT* (LENGTH X)))
  230.     (T (TREEDET X))))
  231.  
  232. (DEFMFUN $IDENT (N) (CONS '($MATRIX) (MXC (DIAGMATRIX N 1 '$IDENT))))
  233.  
  234. (DEFMFUN $DIAGMATRIX (N VAR)
  235.  (CONS '($MATRIX) (MXC (DIAGMATRIX N VAR '$DIAGMATRIX))))
  236.  
  237. (DEFUN DIAGMATRIX (N VAR FN)
  238.        (PROG (I ANS)
  239.          (IF (OR (NOT (EQ (ml-typep N) 'fixnum)) (MINUSP N))
  240.          (IMPROPER-ARG-ERR N FN))
  241.          (SETQ I N)
  242.     LOOP (IF (ZEROP I) (RETURN ANS))
  243.          (SETQ ANS (CONS (ONEN I N VAR 0) ANS) I (f1- I))
  244.          (GO LOOP)))
  245.  
  246. ; ATOMAT GENERATES A MATRIX FROM A MXN ARRAY BY TAKING COLUMNS S TO N
  247.  
  248. (DEFUN ATOMAT (NAME M N S)
  249.   (setq name (get-array-pointer name))
  250.        (PROG (J D ROW MAT)
  251.          (SETQ M (f1+ M) N (f1+ N)) 
  252.     LOOP1(COND ((= M 1) (RETURN MAT)))
  253.          (SETQ M (f1- M) J N)
  254.     LOOP2(COND ((= J S) (SETQ MAT (CONS ROW MAT) ROW NIL) (GO LOOP1)))
  255.          (SETQ J (f1- J))
  256.          (SETQ D (COND (TOP* (MEVAL (LIST (LIST NAME 'array) M J)))
  257.                (T (AREF NAME M J))))
  258.          (SETQ ROW (CONS (OR D '(0 . 1)) ROW))
  259.          (GO LOOP2)))
  260.  
  261. (DEFMFUN $INVERTMX (K) 
  262.   (LET ((*INV* T) *DET* LININD* TOP* MUL* ($RATMX T) (RATMX $RATMX) $RATFAC
  263.     $SPARSE)
  264.     (COND ((ATOM K) ($NOUNIFY '$INVERX) (LIST '(%INVERX) K))
  265.       (T (NEWVARMAT1 (SETQ K (CHECK K)))
  266.          (SETQ K (INVERT1 (REPLIST1 (MCX (CDR K)))))
  267.          (SETQ K (COND ($DETOUT `((MTIMES)
  268.                       ((MEXPT) ,(RDIS (OR *DET* '(1 . 1))) -1)
  269.                       (($MATRIX) ,@(MXC (DISREPLIST1 K)))))
  270.                (T (CONS '($MATRIX) (MXC (DISREPLIST1 K))))))
  271.          (COND ((AND RATMX (NOT $DETOUT))
  272.             (FMAPL1 #'(LAMBDA (X) X) K))
  273.            ((NOT RATMX) ($TOTALDISREP K))
  274.            (T K))))))
  275.  
  276. (DEFUN DIAGINV (AX M)
  277.        (SETQ AX (GET-ARRAY-POINTER AX))
  278.        (COND ($DETOUT (SETQ *DET* 1)
  279.               (DO ((I 1 (f1+ I))) ((> I M))
  280.               (SETQ *DET* (PLCM *DET* (CAR (AREF AX I I)))))
  281.               (SETQ *DET* (CONS *DET* 1))))
  282.        (DO ((I 1 (f1+ I))(ELM))
  283.        ((> I M))
  284.        (SETQ ELM (AREF AX I I))
  285.        (STORE (AREF AX I (f+ M I))
  286.           (COND ($DETOUT (CONS (PTIMES (CDR ELM)
  287.                            (PQUOTIENT (CAR *DET*) (CAR ELM))) 1))
  288.             (T (RATINVERT ELM))))))
  289.  
  290. (DEFUN INVERT1 (K) 
  291.        (PROG (L R G I M N EI* EJ* ONEOFF*) 
  292.          (SETQ L (LENGTH K) I 1) 
  293.          (COND ((= L (LENGTH (CAR K))) NIL)
  294.            (T(MERROR "Non-square matrix in inverse")))
  295.     LOOP (COND ((NULL K) (GO L1))) 
  296.          (SETQ R (CAR K)) 
  297.          (SETQ G (NCONC G (LIST (NCONC R (ONEN I L '(1 . 1) '(0 . 1)))))) 
  298.          (SETQ K (CDR K) I (f1+ I)) 
  299.          (GO LOOP) 
  300.     L1   (SETQ K G)
  301.          (MTOA '*MAT* (SETQ M (LENGTH K)) (SETQ N (LENGTH (CAR K))) K)
  302.          (SETQ K NIL)
  303.          (COND ((DIAGP '*MAT* M) (DIAGINV '*MAT* M)) (T (TFGELI0 '*MAT* M N)))
  304.          (SETQ K (ATOMAT '*MAT* M N (f1+ M)))
  305.          (*REARRAY '*MAT*)
  306.          (RETURN K)))
  307.  
  308. (DEFUN DIAGP (AX M)
  309.    (DECLARE (FIXNUM M ))
  310.    (PROG ((I 0) (J 0))
  311.      (DECLARE (FIXNUM i j))
  312.      (SETQ AX (GET-ARRAY-POINTER AX))
  313.     LOOP1(SETQ I (f1+ I) J 0)
  314.          (COND((> I M) (RETURN T)))
  315.     LOOP2(SETQ J (f1+ J))
  316.          (COND((> J M) (GO LOOP1))
  317.           ((AND (NOT (= I J))(EQUAL (AREF AX I J) '(0 . 1))) NIL)
  318.           ((AND(= I J)(NOT (EQUAL (AREF AX I J) '(0 . 1)))) NIL)
  319.           (T(RETURN NIL)))
  320.      (GO LOOP2)))
  321.  
  322. (DEFUN TFGELI0 (X M N) (COND((OR $SPARSE *DET*) (TFGELI X M N))
  323.                 (T(TFGELI X M N) (DIAGLIZE1 X M N))))
  324.  
  325.     ;  TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE
  326.  
  327. (DEFUN RITEDIV (X M N A)
  328.   (DECLARE(FIXNUM  M N))
  329.   (setq x (get-array-pointer x))
  330.   (PROG ((J 0) (I 0) D ERRRJFFLAG)
  331.     (DECLARE(FIXNUM  i j))
  332.     (SETQ ERRRJFFLAG T)
  333.     (SETQ I M)
  334.   LOOP1 (COND ((ZEROP I) (RETURN NIL)))
  335.     (STORE (AREF X I I) NIL)
  336.     (SETQ J M)
  337.    LOOP (COND ((= J N) (SETQ I (f1- I)) (GO LOOP1)))
  338.     (SETQ J (f1+ J))
  339.     (COND ((EQUAL A 1)
  340.            (STORE (AREF X I J) (CONS (AREF X I J) 1))
  341.            (GO LOOP)))
  342.     (SETQ D (CATCH 'RATERR (PQUOTIENT (AREF X I J) A)))
  343.     (SETQ D (COND (D (CONS D 1)) (T (RATREDUCE (AREF X I J) A))))
  344.     (STORE (AREF X I J) D)
  345.     (GO LOOP)))
  346.  
  347. (DEFUN DIAGLIZE1 (X M N)
  348.   (setq x (get-array-pointer x))
  349.        (PROG NIL
  350.          (COND (*DET* (RETURN (PTIMES *DET* (AREF X M M)))))
  351.          (SETQ *DET* (CONS (AREF X M M) 1))
  352.          (COND ((NOT $DETOUT) (RETURN (RITEDIV X M N (AREF X M M))))
  353.            (T (RETURN (RITEDIV X M N 1))))))
  354.  
  355. ;; Takes an M by N matrix and creates an array containing the elements
  356. ;; of the matrix.  The array is associated "functionally" with the
  357. ;; symbol NAME.
  358. ;; For CL we have put it in the value cell-WFS.  Things still work.
  359.  
  360. (DEFUN MTOA (NAME M N MAT)
  361.        (DECLARE (FIXNUM M N ))
  362.        #+cl
  363.        (proclaim (list 'special name))
  364.        (set name (*array nil t (f1+ M) (f1+ N)))
  365.        (SETQ NAME (get-ARRAY-POINTER NAME))
  366.        (DO ((I 1 (f1+ I))
  367.         (MAT MAT (CDR MAT)))
  368.        ((> I M) NIL)
  369.        (DECLARE(FIXNUM  I))
  370.        (DO ((J 1 (f1+ J))
  371.         (ROW (CAR MAT) (CDR ROW)))
  372.            ((> J N))
  373.            (DECLARE(FIXNUM  J))
  374.            (STORE (AREF NAME I J) (CAR ROW)))))
  375.  
  376.  
  377. (DEFMFUN $ECHELON (X)
  378.   ((LAMBDA ($RATMX) (NEWVARMAT1 (SETQ X (CHECK X)))) T)
  379.   ((LAMBDA (*ECH*)
  380.      (SETQ X (CONS '($MATRIX) (MXC (DISREPLIST1 (ECHELON1 (REPLIST1 (MCX (CDR X)))))))))
  381.    T)
  382.   (COND ($RATMX X) (T ($TOTALDISREP X))))
  383.  
  384. (DEFUN ECHELON1 (X)
  385.   ((LAMBDA (M N)
  386.      (MTOA '*MAT* M N X)
  387.      (SETQ X (CATCH 'RANK (TFGELI '*MAT* M N)))
  388.      (COND ((AND *RANK* X)(THROW 'RNK X))(T (ECHELON2 '*MAT* M N))))
  389.    (LENGTH X) (LENGTH (CAR X))))
  390.  
  391. (DEFUN ECHELON2 (NAME M N)
  392.   (DECLARE (FIXNUM M N ))
  393.   (setq name (symbol-value name))
  394.   (PROG ((J 0) ROW MAT A)
  395.       (DECLARE (FIXNUM J ))
  396.     (SETQ M (f1+ M)) 
  397.      LOOP1(COND ((= M 1) #+MacLisp (*REARRAY NAME) (RETURN MAT)))
  398.     (SETQ M (f1- M) J 0 A NIL)
  399.      LOOP2(COND ((= J N) (SETQ MAT (CONS ROW MAT) ROW NIL) (GO LOOP1)))
  400.     (SETQ J (f1+ J))
  401.     (SETQ ROW (NCONC
  402.             ROW (NCONS
  403.               (COND ((OR(> M J)(EQUAL (AREF NAME M J)  0))
  404.                  '(0 . 1))
  405.                 (A (RATREDUCE (AREF NAME M J)A))
  406.                 (T (SETQ A (AREF NAME M J)) '(1 . 1))))))
  407.     (GO LOOP2)))
  408.  
  409. (DEFUN TRIANG (X)
  410.   ((LAMBDA (M N *TRI*)
  411.      (MTOA '*MAT* M N X) 
  412.      (TFGELI '*MAT* M N)
  413.      (TRIANG2 '*MAT* M N))
  414.    (LENGTH X) (LENGTH (CAR X)) T))
  415.  
  416. (DEFUN TRIANG2 (NAM M N)
  417.   (DECLARE (FIXNUM M N ))
  418.   (SETQ NAM (GET-ARRAY-POINTER NAM))
  419.   (PROG ((J 0) ROW MAT)
  420.     (DECLARE (FIXNUM J))
  421.     (STORE (aref  NAM 0 0) 1)
  422.     (SETQ M (f1+ M)) 
  423.      LOOP1(COND ((= M 1) #+MacLisp (*REARRAY NAM) (RETURN MAT)))
  424.     (SETQ M (f1- M) J 0)
  425.      LOOP2(COND ((= J N) (SETQ MAT (CONS ROW MAT) ROW NIL) (GO LOOP1)))
  426.     (SETQ J (f1+ J))
  427.     (SETQ ROW (NCONC ROW (NCONS
  428.                    (COND ((> M J) '(0 . 1))
  429.                      (T (CONS (AREF NAM M J) 1))))))
  430.     (GO LOOP2)))
  431.  
  432. (DEFMFUN ONEN (N I VAR FILL)
  433.   (PROG (G)
  434.      LOOP (COND ((= I N) (SETQ G (CONS VAR G)))
  435.         ((ZEROP I) (RETURN G)) 
  436.         (T (SETQ G (CONS FILL G))))
  437.     (SETQ I (f1- I))
  438.     (GO LOOP)))
  439.  
  440. (DEFUN TIMEX0 (X Y)
  441.   ((LAMBDA (U V)
  442.      (COND ((AND (NULL U) (NULL V)) (LIST '(MTIMES) X Y))
  443.        ((NULL U) (TIMEX1 X (CONS '($MATRIX) (MCX (CDR V)))))
  444.        ((NULL V) (TIMEX1 Y (CONS '($MATRIX) (MCX (CDR U)))))
  445.        (T (CONS '($MATRIX MULT) (MXC (MULTIPLYMATRICES (MCX (CDR U)) (MCX (CDR V))))))))
  446.    (CHECK1 X) (CHECK1 Y)))
  447.  
  448. (DEFUN TIMEX1 (X Y)
  449.   (SETQ Y (CHECK Y))
  450.   (COND ((NOT $RATMX) (SETQ Y (CDR Y)))
  451.     (T (SETQ X (CDR (RATF X)) Y (REPLIST1 (CDR Y)))))
  452.   (CTIMESX X Y))
  453.  
  454. (DEFUN CTIMESX (X Y)
  455.   (PROG (C)
  456.      LOOP (COND ((NULL Y) 
  457.          (RETURN (CONS '($MATRIX MULT)
  458.                    (MXC (COND ((NOT $RATMX) C) (T (DISREPLIST1 C)))))))) 
  459.     (SETQ C (NCONC C (LIST (TIMESROW X (CAR Y)))) Y (CDR Y))
  460.     (GO LOOP)))
  461.  
  462. (DEFUN MULTIPLYMATRICES (X Y) 
  463.   (COND ((AND (NULL (CDR Y)) (NULL (CDR X)))
  464.      (AND (CDAR X) (SETQ Y (TRANSPOSE Y))))
  465.     ((AND (NULL (CDAR X)) (NULL (CDAR Y)))
  466.      (AND (CDR Y) (SETQ X (TRANSPOSE X)))))
  467.   (COND ((NOT (= (LENGTH (CAR X)) (LENGTH Y)))
  468.      (COND ((AND (NULL (CDR Y)) (= (LENGTH (CAR X)) (LENGTH (CAR Y))))
  469.         (SETQ Y (TRANSPOSE Y)))
  470.            (T (MERROR "incompatible dimensions - cannot multiply")))))
  471.   (COND ((NOT $RATMX) (MULTMAT X Y))
  472.     (T (SETQ X (REPLIST1 X) Y (REPLIST1 Y)) 
  473.        (DISREPLIST1 (MULTMAT X Y))))) 
  474.  
  475. (DEFUN MULTMAT (X Y)
  476.   (PROG (MAT ROW YT ROWX)
  477.     (SETQ YT (TRANSPOSE Y))
  478.      LOOP1(COND ((NULL X) (RETURN MAT)))
  479.     (SETQ ROWX (CAR X) Y YT)
  480.      LOOP2(COND ((NULL Y)
  481.          (SETQ MAT (NCONC MAT (NCONS ROW)) X (CDR X) ROW NIL)
  482.          (GO LOOP1)))
  483.     (SETQ ROW (NCONC ROW (NCONS (MULTL ROWX (CAR Y)))) Y (CDR Y))
  484.     (GO LOOP2)))
  485.  
  486. ;;; This actually takes the inner product of the two vectors.
  487. ;;; I check for the most common cases for speed. '|&*| is a slight
  488. ;;; violation of data abstraction here. The parser should turn "*" into
  489. ;;; MTIMES, well, it may someday, which will break this code. Don't
  490. ;;; hold your breath.
  491.  
  492. (DEFUN MULTL (A B)
  493.        (COND ((EQ $MATRIX_ELEMENT_ADD '|&+|)
  494.           (DO ((ANS (IF (NOT $RATMX) 0 '(0 . 1))
  495.             (COND ((NOT $RATMX)
  496.                    (COND ((EQ $MATRIX_ELEMENT_MULT '|&*|)
  497.                       (ADD ANS (MUL (CAR A) (CAR B))))
  498.                      ((EQ $MATRIX_ELEMENT_MULT '|&.|)
  499.                       (ADD ANS (NCMUL (CAR A) (CAR B))))
  500.                      (T
  501.                       (ADD ANS
  502.                        (MEVAL `((,(GETOPR $MATRIX_ELEMENT_MULT))
  503.                             ((MQUOTE SIMP) ,(CAR A))
  504.                             ((MQUOTE SIMP) ,(CAR B))))))))
  505.                   (T
  506.                    (RATPLUS ANS (RATTIMES (CAR A) (CAR B) T)))))
  507.            (A A (CDR A))
  508.            (B B (CDR B)))
  509.           ((NULL A) ANS)))
  510.          (T
  511.           (MAPPLY (GETOPR $MATRIX_ELEMENT_ADD)
  512.               (MAPCAR #'(LAMBDA (U V)
  513.                  (MEVAL `((,(GETOPR $MATRIX_ELEMENT_MULT))
  514.                       ((MQUOTE SIMP) ,U)
  515.                       ((MQUOTE SIMP) ,V))))
  516.                   A B)
  517.               (GETOPR $MATRIX_ELEMENT_ADD)))))
  518.  
  519. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;       
  520. ;; I leave this for your historical enjoyment. har har.
  521. ;       (PROG (ANS)
  522. ;         (SETQ ANS (COND ((NOT $RATMX) 0) (T '(0 . 1))))
  523. ;    LOOP (COND ((NULL A) (RETURN ANS))) 
  524. ;         (SETQ ANS (COND ((NOT $RATMX)
  525. ;                  (SIMPLUS (LIST '(MPLUS)  ANS  (SIMPTIMES
  526. ;                                 (LIST '(MTIMES)
  527. ;                                   (CAR A)(CAR B))
  528. ;                                 1 T)) 1 T)
  529. ;                  )
  530. ;                 (T (RATPLUS ANS (RATTIMES (CAR A) (CAR B) T)))))
  531. ;         (SETQ A (CDR A) B (CDR B))
  532. ;         (GO LOOP))
  533. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  534.  
  535. (DEFMFUN BBSORT (L FN) (NREVERSE (SORT (copy-top-level L ) FN)))
  536.  
  537. (DEFMFUN POWERX (MAT X) 
  538.        (PROG (N Y) 
  539.          (COND ((NOT (FIXNUMP X))
  540.             (RETURN (LIST '(MNCEXPT SIMP) MAT X)))
  541.            ((= X 1) (RETURN MAT))
  542.            ((MINUSP X)
  543.             (SETQ X (MINUS X) MAT ($INVERTMX MAT))
  544.             (COND ($DETOUT
  545.                (RETURN (LET ((*INV* '$DETOUT))
  546.                     (MUL2*
  547.                      (POWER* (CADR MAT) X)
  548.                      (FMAPL1 #'(LAMBDA (X) X)
  549.                          (POWERX (CADDR MAT) X)))))))))
  550.          (NEWVARMAT1 (SETQ MAT (CHECK MAT)))
  551.          (SETQ N 1 MAT (MCX (CDR MAT)) Y MAT) 
  552.     LOOP (IF (= N X)
  553.          (LET (($SCALARMATRIXP (IF (EQ $SCALARMATRIXP '$ALL) '$ALL)))
  554.            (RETURN (SIMPLIFY (CONS '($MATRIX MULT) (MXC Y))))))
  555.          (SETQ Y (MULTIPLYMATRICES Y MAT) N (f1+ N)) 
  556.          (GO LOOP))) 
  557.  
  558. ;; The following $ALGEBRAIC code is so that 
  559. ;; RANK(MATRIX([1-SQRT(5),2],[-2,1+SQRT(5)])); will give 1.
  560. ;; - JPG and BMT
  561.  
  562. (DEFMFUN $RANK (X)
  563.   (LET ((*RANK* T) ($RATMX T) ($ALGEBRAIC $ALGEBRAIC))
  564.     (NEWVARMAT1 (SETQ X (CHECK X)))
  565.     (AND (NOT $ALGEBRAIC) (ORMAPC #'ALGP VARLIST) (SETQ $ALGEBRAIC T))
  566.     (SETQ X (REPLIST1 (MCX (CDR X))))
  567.     (MTOA '*MAT* (LENGTH X) (LENGTH (CAR X)) X)
  568.     (TFGELI '*MAT* (LENGTH X) (LENGTH (CAR X)))))
  569.  
  570. (DEFUN REPLACEROW (I Y X)
  571.  (IF (= I 1)
  572.      (cons y (cdr x)) ;(NCONC (LIST Y) (CDR X))
  573.      (cons (car x) (REPLACEROW (f1- I) Y (CDR X)))
  574.      ;(NCONC (LIST (CAR X)) (REPLACEROW (f1- I) Y (CDR X)))
  575.      ))
  576.  
  577. (DEFUN TIMESROW (Y ROW)
  578.  (PROG (ANS)
  579.        (COND ((AND $RATMX (ATOM Y) Y) (SETQ Y (CDR (RATF Y)))))
  580.   LOOP (COND ((NULL ROW) (RETURN ANS)))
  581.        (SETQ ANS (NCONC ANS (LIST (COND ((NOT $RATMX)
  582.                      (SIMPTIMES
  583.                       (LIST '(MTIMES) Y (CAR ROW)) 1 NIL))
  584.                     (T (RATTIMES Y (CAR ROW) T))))))
  585.        (SETQ ROW (CDR ROW))
  586.        (GO LOOP)))
  587.  
  588. (DEFMFUN $TRIANGULARIZE (X) 
  589.   ((LAMBDA ($RATMX) (NEWVARMAT1 (SETQ X (CHECK X)))) T)
  590.   (SETQ X (CONS '($MATRIX) (MXC (DISREPLIST1 (TRIANG (REPLIST1 (MCX (CDR X)))))))) 
  591.   (COND ($RATMX X) (T ($TOTALDISREP X))))
  592.  
  593. (DEFMFUN $COL (MAT N)
  594.   (CONS '($MATRIX) (MXC (TRANSPOSE (LIST (NTHCOL (MCX (CDR (CHECK MAT))) N)))))) 
  595.  
  596. (DEFUN DELETECOL (N X)
  597.        (PROG (M G)
  598.          (SETQ M X)
  599.     LOOP (COND ((NULL M) (RETURN G)))
  600.          (SETQ G (NCONC G (NCONS (DELETEROW N (CAR M)))) M (CDR M))
  601.          (GO LOOP)))
  602.  
  603. (DEFUN DELETEROW (I M) 
  604.   (COND ((OR (NULL M) (LESSP I 0)) (MERROR "Incorrect index - MATRIX"))
  605.     ((= I 1) (CDR M)) 
  606.     (T (CONS (CAR M) (DELETEROW (f1- I) (CDR M)))))) 
  607.  
  608. (DEFMFUN $MINOR (MAT M N) (CONS '($MATRIX) (MXC (MINOR M N (MCX (CDR (CHECK MAT)))))))
  609.  
  610. (DEFUN MINOR (I J M) (DELETECOL J (DELETEROW I M))) 
  611.  
  612. (DEFMFUN $ROW (MAT M) (CONS '($MATRIX) (MXC (LIST (ITH (MCX (CDR (CHECK MAT))) M)))))
  613.  
  614. (DEFMFUN $SETELMX (ELM M N MAT) 
  615.   (COND ((NOT (AND (INTEGERP M) (INTEGERP N) ($MATRIXP MAT)))
  616.      (MERROR "Wrong arg to SETELMX"))
  617.     ((NOT (AND (> M 0) (> N 0) (> (LENGTH MAT) M) (> (LENGTH (CADR MAT)) N)))
  618.      (MERROR "No such entry - SETELMX")))
  619.   (RPLACA (NCDR (CAR (NCDR MAT (f1+ M))) (f1+ N)) ELM) MAT) 
  620.  
  621. ;;; Here the function transpose can actually do simplification of
  622. ;;; its argument. TRANSPOSE(TRANSPOSE(FOO)) => FOO.
  623. ;;; If you think this is a hack, well, realize that the hack is
  624. ;;; actually the fact that TRANSPOSE can return a noun form.
  625.  
  626. (DEFMFUN $TRANSPOSE (MAT)
  627.   (COND ((NOT (MXORLISTP MAT))
  628.      (COND ((AND (NOT (ATOM MAT))
  629.              (EQ (CAAR MAT) '%TRANSPOSE))
  630.         (CADR MAT))
  631.            (($SCALARP MAT) MAT)
  632.            ((MPLUSP MAT)
  633.         `((MPLUS) .,(MAPCAR #'$TRANSPOSE (CDR MAT))))
  634.            ((MTIMESP MAT)
  635.         `((MTIMES) .,(MAPCAR #'$TRANSPOSE (CDR MAT))))
  636.            ((MNCTIMESP MAT)
  637.         `((MNCTIMES) .,(NREVERSE (MAPCAR #'$TRANSPOSE (CDR MAT)))))
  638.            ((MNCEXPTP MAT)
  639.         (LET (((MAT POW) (CDR MAT)))
  640.           `((MNCEXPT) ,($TRANSPOSE MAT) ,POW)))
  641.            
  642.            (T ($NOUNIFY '$TRANSPOSE) (LIST '(%TRANSPOSE) MAT))))
  643.     (T
  644.      (LET ((ANS (TRANSPOSE (MCX (CDR (CHECK MAT))))))
  645.        (COND ($MATRIX_ELEMENT_TRANSPOSE
  646.           (SETQ ANS (MAPCAR #'(LAMBDA (U)
  647.                     (MAPCAR #'TRANSPOSE-ELS
  648.                         U))
  649.                     ANS))))
  650.        `(($MATRIX) . ,(MXC ANS))))))
  651.  
  652. ;;; THIS IS FOR TRANSPOSING THE ELEMENTS OF A MATRIX
  653. ;;; A hack for Block matricies and tensors.
  654.  
  655. (DEFUN TRANSPOSE-ELS (ELEM)
  656.        (COND ((EQ $MATRIX_ELEMENT_TRANSPOSE '$TRANSPOSE)
  657.           ($TRANSPOSE ELEM))
  658.          ((EQ $MATRIX_ELEMENT_TRANSPOSE '$NONSCALARS)
  659.           (COND (($NONSCALARP ELEM)
  660.              ($TRANSPOSE ELEM))
  661.             (T ELEM)))
  662.          (T
  663.           (MEVAL `((,(GETOPR $MATRIX_ELEMENT_TRANSPOSE)) ((MQUOTE SIMP) ,ELEM))))))
  664.  
  665.  
  666. (DEFMFUN $SUBMATRIX NARGS
  667.  (PROG (R C X)
  668.        (SETQ X (LISTIFY NARGS))
  669.   L1   (COND ((NUMBERP (CAR X)) (SETQ R (CONS (CAR X) R) X (CDR X)) (GO L1)))
  670.        (SETQ C (NREVERSE (BBSORT (CDR X) '>)) R (NREVERSE (BBSORT R '>)))
  671.        (SETQ X (MCX (CDAR X)))
  672.   L2   (COND ((NULL R) (GO B)) (T (SETQ X (DELETEROW (CAR R) X))))
  673.        (SETQ R (CDR R))
  674.        (GO L2)
  675.   B    (COND ((NULL C) (RETURN (CONS '($MATRIX) (MXC X)))))
  676.        (SETQ X (DELETECOL (CAR C) X) C (CDR C))
  677.        (GO B)))
  678.  
  679.  
  680. (defun $LIST_MATRIX_ENTRIES (m)
  681.   (or ($matrixp m) (error "not a matrix"))
  682.   (cons (if (null (cdr m)) '(mlist) (caadr m))
  683.     (sloop for row in (cdr m) append (cdr row))))
  684.  
  685. ; Undeclarations for the file:
  686. #-NIL
  687. (DECLARE-TOP(NOTYPE NN LEN))
  688.